% writen by ruiguo
%i_SSize was meaned two images have same  size

function [ ravg, rstd] = readavgsdofphantom( file_path,Ref_mask_map,t1_t2,i_SSize)

load( file_path );
slice_num = 2;
if( size(T1,3) < 2 )
    slice_num = 1 ;
end
Tx_map = squeeze( T1(:,:,slice_num,t1_t2));
Tx_map(isnan(Tx_map)) = 0;
%Tx_map = fliplr(Tx_map);
Bineary_map = zeros( size(Tx_map));
%for weight
Tx_map(Tx_map<((t1_t2==1)*100+(t1_t2==2)*10)) = 0;
Bineary_map(Tx_map>0) = 1;

Label_Bineary_map = zeros( size(Tx_map));
[Label_Bineary_map, object_num] = bwlabel(Bineary_map,4);
Theta_arr = linspace(0 , 2*pi , 360);
Mask_num = size( Ref_mask_map, 3 );
Average_arr = zeros( 1, Mask_num);
Std_arr = zeros( 1, Mask_num);
if(isempty(i_SSize))
    
    
    for object_ix = 1: object_num
        Temp_map =  zeros( size(Tx_map));
        Temp_map(Label_Bineary_map==object_ix) = 1;
        [row_iy, col_ix,value]=find(Temp_map>0);
        area_size = sum( sum( Temp_map ));
        if( area_size < 50  )
            Label_Bineary_map(Label_Bineary_map==object_ix) = 0;
        end
    end
    Label_Bineary_map(Label_Bineary_map>0) = 1;
    [Label_Bineary_map, object_num] = bwlabel(Label_Bineary_map,4);
    for object_ix = 1: object_num
        Temp_map =  zeros( size(Tx_map));
        Temp_map(Label_Bineary_map==object_ix) = 1;
        [row_iy, col_ix,value]=find(Temp_map>0);
        area_size = sum( sum( Temp_map ));
        Center_iy = mean( row_iy );
        Center_ix = mean( col_ix );
        Obj_num_ix = 0;
        Mask_center_dis = zeros(1,Mask_num);
        for mask_ix  = 1 : Mask_num
            mask_map = Ref_mask_map(:,:,mask_ix);
            [mask_row_iy, mask_col_ix,value]=find(mask_map>0);
            mask_Center_iy = mean( mask_row_iy );
            mask_Center_ix = mean( mask_col_ix );
            Mask_center_dis( 1, mask_ix ) = sqrt( ( mask_Center_iy - Center_iy  ).^2+( mask_Center_ix - Center_ix ).^2 );
        end
        [min_row, min_col, value] =  find( Mask_center_dis == min(Mask_center_dis));
        if(min_col(1,1) > 0 )
            Obj_num_ix = min_col(1,1) ;
            if(size(Tx_map,1)~=size(Ref_mask_map,1 ))
                
                Row_len = max( row_iy  ) - min( row_iy );
                Col_len = max( col_ix ) - min( col_ix );
                R = min( [Row_len Col_len ])*0.40;     %RuiGuo: Changed it from 0.45 to 0.25 on 2017-07-15
                Y_arr = R*sin(Theta_arr)+ Center_iy ;
                X_arr = R*cos(Theta_arr) + Center_ix ;
                %mesh grid
                [mesh_x,mesh_y] = meshgrid([1:size(Tx_map,2)],[1:size(Tx_map,1)]);
                obj_area = inpolygon(mesh_x,mesh_y,X_arr,Y_arr);
                %             [plant_i,plant_j] = find(obj_area);
                %             obj_area = zeros(size(T1_mapping),size(T1_mapping));
            else
                obj_area = squeeze( Ref_mask_map(:,:,Obj_num_ix));
                obj_area = obj_area.*Tx_map ;
            end
            Average_arr(1,Obj_num_ix) = mean( obj_area(obj_area >0) );
            Std_arr( 1, Obj_num_ix )= std( obj_area(obj_area >0) );
            
            obj_area( obj_area >3000) = 0;
            Average_arr(1,Obj_num_ix) =round(10*mean( obj_area(obj_area >0)) )/10;
            Std_arr( 1, Obj_num_ix )=round(10*std( obj_area(obj_area >0)))/10 ;
        end
    end
else

    object_num = Mask_num;
    for object_ix = 1: object_num
        Obj_num_ix = object_ix;
        obj_area = squeeze( Ref_mask_map(:,:,Obj_num_ix));
        if ~isequal( size(obj_area ), size(Tx_map) )
            obj_area = imresize(obj_area, [size(Tx_map,1) size(Tx_map,2)], 'nearest');        
        end

        if(0)
            ROI_temp = edge(sum(Ref_mask_map,3));
            ROI_temp = ROI_temp.*((t1_t2==1)*(-3000)+(t1_t2==2)*(-500));
            ROI_temp = Tx_map+ROI_temp;
            f = figure('Name','ROI of refference', 'visible','on');
            imagesc(ROI_temp,(t1_t2==1)*[ 0 2000]+(t1_t2==2)*[0 300]),colormap(jet);
            axis normal;
            axis off;
        end

        obj_area = obj_area.*Tx_map ;
        
        Average_arr(1,Obj_num_ix) = mean( obj_area(obj_area >0) );
        Std_arr( 1, Obj_num_ix )= std( obj_area(obj_area >0) );
        
       obj_area( obj_area > 3000 ) = 0;
        Average_arr(1,Obj_num_ix) =round(10*mean( obj_area(obj_area >0)) )/10;
        Std_arr( 1, Obj_num_ix )=round(10*std( obj_area(obj_area >0)))/10 ;
    end
end
ravg = Average_arr;
rstd = Std_arr ;
end